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Abstract 

After developing an appropriate iteration procedure for the determination 
of the parameters, the method of simulated tempering has been successfully 
applied to the 2D Ising spin glass. The reduction of the slowing down is 
comparable to that of the multicanonical algorithm. Simulated tempering 
has, however, the advantages to allow full vectorization of the programs and 
to provide the canonical ensemble directly. 
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I. Introduction 



A better theoretical understanding of spin glasses is still highly desirable. To 
make progress within this respect more efficient simulation algorithms are needed. 
In particular at low temperatures conventional simulations suffer from severe slow- 
ing down due to energy barriers. Recently Berg and Celik [|lj have been able to 
reduce the slowing down considerably by applying the multicanonical method [Q. 
Shortcomings of this method are that the computer programs cannot be vectorized 
and that additional calculations (with complicated error estimates) are needed to 
obtain the canonical ensemble. 

An alternative is offered by the method of simulated tempering which has been 
introduced by Marinari and Parisi || for the random-field Ising model. It works by 
associating a set of inverse temperatures to an additional dynamical variable of the 
Monte Carlo scheme. Because then the system is warmed up from time to time one 
gets a chance to avoid the effective ergodicity breaking M observed in conventional 
approaches. By this method at each of the temperatures one gets the canonical 
ensemble directly and there is nothing that prevents vectorization of the programs. 

In these simulations involving the joint probability distribution of the spins and 
of the additional variable one has the freedom to choose the set of temperatures 
and the probability distribution of the additional variable. To fix these quantities 
in an optimal way is crucial for the efficiency of the method. In the first applica- 
tion to the random-field Ising model on small lattices this appeared to be not very 
demanding [[|. However, for spin glasses, in particular for larger lattices and low 
temperatures, we find that the appropriate determination of these parameters is far 
from straightforward. 

In the present paper we discuss the related issues in detail and develop a systematic 
procedure which allows to fix the respective parameters in an appropriate way. We 
apply our method to the 2D Ising spin glass and show that an efficient algorithm is 
obtained. 

We use the Edwards-Anderson Hamiltonian with two replicas, 

H(s,t) = -^2 Jij(siSj + Utj) , (1.1) 

where the sum goes over nearest neighbors, and Sj = ±1, £, = ±1. Thus we are able 
to measure the order parameter [|, [5], |J 

g = X> f ti • (1.2) 

i 

We investigate samples for which is put equal to +1 and —1 with equal proba- 
bility. 
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To each value of the additional dynamical variable n of simulated tempering, which 
takes the values 1, 2, . . . , N, an inverse temperature (3 n with f3\ < fii < . . . < (3n is 
associated. The partition function of the joint system then becomes 

Z =J2^v(-PnH(s,t)+g n ) . (1.3) 

n,s,t 

The parameters g n reflect the freedom in the choice of the probability distribution 
p(n) of the variable n. 

Two steps are necessary in the calculations. In the first step the parameters g n and 
P n are determined. Then in the second step, using these parameters, the simulations 
are performed and the physical observables are measured. The dynamical variables 
ti, n are all updated by standard Metropolis || techniques. Slowing down is 
measured in terms of the ergodicity time te, which is defined as the mean time the 
system needs to travel from n = 1 to iV and back. 

In Sec. II the iteration procedures for the determination of the parameters is 
described. Sec. Ill contains a discussion of the properties of this procedure. In 
Sec. IV some numerical results are presented. 



II. Iteration procedure 

To choose p(n) constant as proposed in Ref. ||, i.e. to visit the states n with equal 
probability, appears reasonable. In a similar approach 0, in which the additional 
dynamical variables is the number of components in the Potts model, it has been 
tested that such a choice is optimal. Constant p(n) in terms of the g n means that 

g n = -\nZ(P n ) (2.1) 

with the spin-glass partition function 

Z(/3) = £exp(-/3#(s,*)) . (2.2) 

s,t 

It should be noted that g n depends on (3 n only. 

Our choice here is to come sufficiently close to (|2.1|) . Because of the exponential 
dependence of p(n) on g n the deviations from ( |2.1p must be small in order that all 
states n are appropriately reached. We use two methods to calculate the g n . The 
first one is to replace g n by g n — ln(Np(n)) in subsequent iterations such that one 
ultimately arrives at constant p(n). The second method is to apply reweighting || 
of the obtained distribution to get the ratio of the Z at neighboring n and thus by 
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( |2.1| ) the differences of the respective new g n . The new g n then are calculated from 
these differences. 

To get hints on a reasonable choice of the (3 n one observes that the mean value of 
the logarithm of the acceptance probability exp(— Af3H + Ag) to order (A/?) 3 is 

H 2 ) - (H) 2 ) (A/?) 2 . (2.3) 

On this basis one may require H to choose f3 n in such a way that (|2.3| ) gets as 
constant as possible. Then, due to (H 2 ) — (H) 2 ~ L 2 at fixed f3, it follows that 
the number N of temperature values should scale like iV ~ L on the L x L lattice. 
Further, this indicates that at larger (3, where (H 2 ) — (H) 2 is small, the distances 
between the /^-values should increase. Unfortunately at low temperatures (large (3) 
the requirement to make ( |2.3|) constant cannot be used to determine Aj3 because 
the measurements are likely to arrive at (H 2 ) — (H) 2 = 0, which would result in a 
breakdown of an iterative determination of the f3 n . 

We propose to calculate the (3 n by using a map which has an attractive fixed point. 
The map which we have constructed and applied has the property that if the fixed 
point is reached in a sequence of mappings the effective stay time at a state n 
gets constant. This quantity is defined by 

eff f 2 Tn ' n 1> (9 A] 

n ~ \r n , 2 < n < iV - 1 1 ' 

where the stay time r n is the mean time which elapses between entering state n and 
leaving it. To perform the map we first compute the auxiliary quantities 

a n = (P n+ i-Pn)/K{{ + < ff ), n = l,...,N-l , (2.5) 

N-l 

c = J2a n . (2.6) 

n=l 

Then we get new values /3' n by 

PI = Pi , (2.7) 

P' n = t^^ /""^ 1 , n = 2,...,N . (2.8) 

c 

Since /3i and [3n are mapped to themselves the region of (3 covered does not change. 
The fixed point gets attractive because a large stay time implies a smaller (3- 
difference and thus a smaller stay time in the next iteration step. The numerical 
results indicate that by the sequence of mappings also ( |2.3|) gets very close to a 
constant value. 

Our procedure, which determines the parameters in an iterative way, is 



started with equidistant (3 n and with g n estimated by the relation [|TTJ g n 
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— L 2 (a/3 n + 6exp(— k(3)) + c. Each iteration involves four steps: 1) An appropriate 
number (typically some multiple of the expected ergodicity time) of Monte Carlo 
sweeps are performed to obtain the necessary data, as e.g. p(n) and r n . 2) New g n 
are computed by one of the two methodes indicated above. 3) The map described 
above is used to obtain new /3 n . 4) New g n which are related to these /3 n are de- 
termined by spline interpolation [|l(J (this step is dictated by the fact that g n is a 
function of (3 n ). Criteria for the termination of the iteration procedure are discussed 
in Sec. III. 



III. Properties of the procedure 

To study the impact of the number N of /3-values, we have measured the ergod- 
icity time as a function of N. Fig. 1 shows the respective results for lattices 12 x 12 
and 24 x 24. These results suggest to avoid the steep increase at small N by choosing 
N slightly above the minimum. Thus we put N = 1.25L in our simulations. This is 
in accordance with the remark in Sec. II that one should have N ~ L and it turns 
out to be appropriate on all lattices considered. 

On smaller lattices we have used the first method to get the new f3 n described 
in Sec. II. On large lattices the second one involving reweighting, which performs 
better for larger N, has been applied. 

For lattices with L > 24 we find effects caused by the occurrence of of metastable 
spin-glass states. Since g n = /3 n f n , where f n is the free energy, one should have 

^ = H for /5»1 , (3.1) 

where H is the ground state energy. In our calculations (g^ — 9n-i)/(Pn ~ Pn-i) 
turns out to be, in fact, with good accuracy an integer number. However, one 
actually relies on the lowest energy which has been reached in the Monte Carlo 
sweeps, which might be larger than H and thus lead to a wrong determination of 
the parameters. 

To see which serious trouble this can cause suppose that 

f "f- 1 =H + AH (3.2) 

PN ~ PN-l 

with AH > (where AH can take the values 4,8,...). The probability to get 
from state iV (related to the largest (3) to state N — 1 is min(p^, 1) with = 
exp (— (/3jv-i — Pn)H + (gN-i — 9n)) which after inserting (|3.2|) becomes 



p A = ex P ((p N -p N _ 1 )(H-H -AH)) . (3.3) 
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From ( p.3| ) it is seen that if in the course of the simulations a true ground state with 
H = Ho is reached, the probability for the transition from N to N — 1 becomes 
extremely small, i.e. the system gets trapped at low temperature. 



Therefore special care is needed to avoid a wrong determination of type ( |3.2|) of 
the parameters. For this reason we have done a large number of iterations in each of 
which we have determined the set of parameters and an estimate of the ergodicity 
time te- Out of these sets we have selected the set of parameters associated to the 
smallest te which in addition has satisfied the requirement 



9n — 9n-i tj 



Pn — Pn-i 



<e . (3.4) 



In ( |3.4| ) for H the lowest energy reached in all previous iterations has been inserted, 
which has turned out to be a workable concept. Because, as already pointed out, 
for Ag/A/3 we find integers with good accuracy, e = 0.1 is a reasonable bound. The 
iterations have been terminated if after an appropriate time no set of parameters 
with still lower te and respecting (|3.2|) has occurred. 

On the largest lattices which we have considered to reach sufficient accuracy of the 
parameters in the iterations has turned out to be cumbersome. If the parameters 
used in a simulation are not accurate enough, the sweeps happen to get restricted 
to a subinterval of the n-interval which spoils the calculation. 



IV. Numerical results 

Our simulations have been performed on 2D lattices with periodic boundary 
conditions of sizes L — 4, 12, 24, 32, 48 for ten samples of Jjj-configurations in each 
case. This has been done for 0.3 < (3 < 3.5 putting N = 1.25L. The ground-state 
energy density 

Eo = j2 m i n ^~ ^ J ij s i s j\ , (4- 1 ) 

the distribution P(q) of the order parameter (|1.2|), and the ergodicity time te have 
been determined. From the obtained data the moments (g 2 ) and (q 4 ) of P(q) and 
the Binder parameter 

have been calculated. The results of this for (3 = 3.5 are presented in Tables 1-5. 
For the samples the errors of the listed quantities are statistical ones including the 
effect of autocorelations (the Eq are exact). The errors of the mean values are the 
ones derived from the fluctuations of the sample values. 



^=n(3-^] ■ (4-2) 
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Our results for the mean value of the ground-state energy within errors agree with 
the numbers of other authors |], 11], 12| . The results for (q 2 ) and B q are compatible 



with those given by Berg and Celik jrj, though they are generally slightly larger. A 
possible reason for this is that we have measured at a lower temperature than they 
did. From our tables it is seen that there are strong dependences on the particular 
sample as is expected because of the lack of self-averaging ||, ||, || . 

Figs. 2 and 3 show typical results for the distribution P(q). The invariance of the 
Hamiltonian under the transformation Sj — > Sj, i; — > —t% requires P(q) = P{—q) as 
is observed. The scaling law [ 13 1 P (q) = L 01 P{qL - 1 ) with an universal function P 



was verified (due to the low number of samples with rather large errors). 

The dependence of the ergodicity time on the lattice size L is depicted in Fig. 4 
and compared to the multicanonical data given by Berg and Celik We get the 
dependence 

t e ~ £ 4 - 27 ( 8 ) . (4.3) 

A fit of the form te ~ exp(A;L) is not possible. Our ergodicity times are comparable 
with those of Ref. (where the dynamical criticial exponent 4.4(3)) is obtained). 
However, because our computer programs can be fully vectorized in terms of CPU 
times we gain a large factor. 
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Table I 



Numerical results for 4x4 lattice. 



# 




E 


(<7 2 > 


<<7 4 > 


B q 


1 


70.53(8) 


-1.125 


0.2419(6) 


0.1247(5) 


0.435(9) 


2 


92.00(11) 


-1.500 


0.88256(29) 


0.7927(5) 


0.9912(7) 


3 


99.84 13 


-1.375 


0.7986 4) 


0.6642(7) 


0.9793(10) 


4 


78.39(9) 


-1.250 


0.4477(6) 


0.2767(6) 


0.8099(35) 


5 


64.05(7) 


-1.000 


0.1746(4) 


0.07206(30) 


0.318(11) 


6 


86.01(10) 


-1.750 


1.00000(1) 


1.00000(1) 


1.0000(1) 


7 


85.88(10) 


-1.375 


0.6620(6) 


0.4923(7) 


0.9383(18) 


8 


212.2(4) 


-1.250 


0.99983(4) 


0.99979(5) 


1.0000(1) 


9 


89.02(11) 


-1.375 


0.6956(4) 


0.5174(6) 


0.9653(13) 


10 


100.28(13) 


-1.250 


0.5237(10) 


0.4372(10) 


0.703(5) 


Mean 


97(13) 


-1.33(7) 


0.64(9) 


0.54(10) 


0.81(8) 



Table II 



Numerical results for 12 x 12 lattice. 



# 




E 


(? 2 ) 


(? 4 ) 


B q 


1 


1689(14) 


-1.3611 


0.351(4) 


0.1660(28) 


0.826(15) 


2 


4576(61) 


-1.3611 


0.458(5) 


0.248(5) 


0.911(25) 


3 


2818(30) 


-1.4444 


0.778(2) 


0.610(2) 


0.996(4) 


4 


3509(42) 


-1.3611 


0.481(7) 


0.311(7) 


0.83(4) 


5 


8421(153) 


-1.3333 


0.592(4) 


0.360(4) 


0.986(12) 


6 


1219(8) 


-1.3750 


0.2250(27) 


0.0872(17) 


0.64(4) 


7 


3944(49) 


-1.3472 


0.5097(32) 


0.2793(33) 


0.962(13) 


8 


1271(9) 


-1.3472 


0.2732(29) 


0.1210(19) 


0.690(30) 


9 


2722(29) 


-1.4861 


0.7842(31) 


0.640(4) 


0.980(8) 


10 


1360(10) 


-1.4444 


0.302(4) 


0.1656(31) 


0.59(4) 


Mean 


3153(695) 


-1.386(17) 


0.48(6) 


0.30(6) 


0.84(5) 
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Table III 



Numerical results for 24 x 24 lattice. 



# 


TE 


E 


(<7 2 > 


(? 4 ) 


B q 


1 


31963(516) 


-1.3958 


0.472(5) 


0.255(5) 


0.929(24) 


2 


55208(1175) 


-1.4271 


0.430(9) 


0.225(7) 


0.89(5) 


3 


184218(7176) 


-1.3993 


0.532(38) 


0.314(33) 


0.94(14) 


4 


76188(1905) 


-1.4028 


0.576(3) 


0.342(4) 


0.985(11) 


5 


108775(3289) 


-1.4167 


0.619(2) 


0.393(2) 


0.986(7) 


6 


43150(826) 


-1.3958 


0.329(5) 


0.130(3) 


0.898(30) 


7 


97243(2810) 


-1.4271 


0.524(5) 


0.304(5) 


0.946(18) 


8 


192333(7852) 


-1.3958 


0.401(31) 


0.250(21) 


0.72(17) 


9 


109485(3291) 


-1.4063 


0.599(6) 


0.382(7) 


0.968(20) 


10 


40648(751) 


-1.3889 


0.328(8) 


0.167(6) 


0.72(7) 


Mean 


93921(18067) 


-1.406(4) 


0.481(34) 


0.276(28) 


0.900(31) 



Table IV 



Numerical results for 32 x 32 lattice. 



# 


10- 5 T E 


E 


<<7 2 > 


(? 4 ) 


B q 


1 


3.73(10) 


-1.4082 


0.500(15) 


0.271(14) 


0.96(6) 


2 


15.0(6) 


-1.3906 


0.43(4) 


0.208(33) 


0.93(20) 


3 


11.0(5) 


-1.3867 


0.348(13) 


0.177(7) 


0.77(8) 


4 


1.147(16) 


-1.4277 


0.229(5) 


0.099(3) 


0.55(7) 


5 


2.73(5) 


-1.4121 


0.443(6) 


0.217(5) 


0.949(28) 


6 


8.74(34) 


-1.4082 


0.447(20) 


0.250(14) 


0.88(9) 


7 


5.79(17) 


-1.3945 


0.658(6) 


0.438(8) 


0.994(10) 


8 


3.45(8) 


-1.4063 


0.162(14) 


0.049(7) 


0.57(29) 


9 


5.83(16) 


-1.3965 


0.576(11) 


0.343(10) 


0.983(34) 


10 


2.60(5) 


-1.4121 


0.490(6) 


0.268(5) 


0.941(23) 


Mean 


6.0(14) 


-1.404(4) 


0.43(5) 


0.232(35) 


0.85(5) 
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Table V 



Numerical results for 48 x 48 lattice. 





io-% 


Eo 


(q 2 ) 


(<7 4 > 


B q 


1 


4.83(27) 


-1.3967 


0.304(8) 


0.117(5) 


0.87(6) 


2 


12.0(10) 


-1.3967 


0.054(10) 


0.004(3) 


0.74(79) 


3 


5.5(3) 


-1.4115 


0.35(1) 


0.131(7) 


0.98(6) 


4 


2.81(15) 


-1.4054 


0.143(12) 


0.021(6) 


0.99(24) 


5 


2.06(11) 


-1.3906 


0.23(5) 


0.086(26) 


0.71 57 


6 


6.3(4) 


-1.4063 


0.48(22) 


0.26(19) 


0.93(92) 


7 


6.02(34) 


-1.4019 


0.549(23) 


0.316(20) 


0.98(8) 


8 


1.00(4) 


-1.3845 


0.370(7) 


0.146(5) 


0.97(4) 


9 


9.2(7) 


-1.4115 


0.42(4) 


0.18(4) 


0.97(20) 


10 


5.1(3) 


-1.3976 


0.28(4) 


0.085(27) 


0.94(35) 


Mean 


5.5(10) 


-1.4003(28) 


0.32(4) 


0.135(31) 


0.91(3) 
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Figure Captions 



Fig. 1. Ergodicity time r E as a function of N for lattices 

12 x 12 (dots) and 24 x 24 (crosses). 
Fig. 2. Order-parameter distribution P{q,j3) for sample 5 

on the 32 x 32 lattice. 
Fig. 3. Order-parameter distribution P(q,(3) for sample 3 

on the 32 x 32 lattice. 
Fig. 4. Dependence of te on the lattice size L for simulated 

tempering (dots, fit solid line) and the multicanonical method 

(crosses, fit dashed line). 
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